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Abstract 

We characterize numerically the properties of the phase transition of the three 
dimensional Ising spin glass with Gaussian couplings and of the low temperature 
phase. We compute critical exponents on large lattices. We study in detail the overlap 
probability distribution and the equilibrium overlap-overlap correlation functions. 
We find a clear agreement with off-equilibrium results from previous work. These 
results strongly support the existence of a continuous spontaneous replica symmetry 
breaking in three dimensional spin glasses. 



PACS mumbers: 75.50.Lk, 05.50.+q, 64.60.Cn 



1 Introduction 



Recent numerical simulations [|T], § have given a strong numerical evidence of the existence 
of a spin glass phase transition in the 3D Edwards Anderson spin glass (for a critical point 
of view see || f|, [5]]). The first studies of such models are today 15 years old (for a recent 
review see ||), and the Replica Symmetry Breaking (RSB) mean field solution |7j is found 
to describe the most of the properties observed in finite dimensional models. The first issue 
is the study of the probability distribution of the order parameter P(q), but the innovative 
features of the RSB solution inspire many questions that can be answered numerically. 

Recently many of these issues have indeed been analyzed numerically. Correlation 
functions and block overlaps of the 3D spin glass have been shown in M to have a mean 
field like behavior, the 4D model has been studied in detail in fllQl , and the role of the 



upper critical dimension, D^ = 6, has been analyzed in [0]|. It is also remarkable that 
a deep relation among static and dynamic behavior has be shown to be valid, for the 3D 
and 4D models, in At last it is worth to remind the reader that results about the 



transition in field in finite dimensional models have been obtained in [|13|, [14| . 

Here we try to push these results to the limit allowed by the computers we can use 
today, and by the numerical improved algorithms we use for our simulations. We consider 
the 3D spin glass with gaussian couplings, that make safer the approach to the low T region 
(because of the absence of accidental degeneracy) and allow to check for universality when 
comparing to the results of 0, [|. After defining our model and our statistical quantities 
and giving details about the numerical simulations we determine critical exponents by 
analyzing the overlap Binder cumulant and the spin glass susceptibility in section (|j) . We 
analyze in detail the behavior of P(q) in (|3p, by checking in detail the presence of a mean 
field like behavior. We analyze the behavior of a simple function in (|(|). In section (^) we 
are able to compute equilibrium correlation function, for all states and for the zero overlap 
sector of the theory, and to show that these results coincide with the one obtained by using 
an extrapolation in |J. 



2 The Model 

We consider a three dimensional {3D) Edwards Anderson spin glass model with Gaussian 
quenched random couplings J. The model is defined on a simple cubic lattice with periodic 
boundary conditions. The Hamiltonian of the system is 

H = - a iJijVj , (1) 

<ij> 

where by < ij > we indicate that the sum runs over couples of first neighboring sites. The 
Jij are quenched Gaussian variables with zero mean and unit variance. We consider two 
real replicas of the system (in the same realization of the disorder). We define the overlap 
among the two real replicas a and 7 at site i as 

tig = °>1 , (2) 

(where the subscript J reminds us that we are in a fixed realization of the quenched disorder 
J), and a total overlap 



2 



(3) 



where we will frequently ignore the superscripts, and denote it by q. Its probability distri- 
bution for a given sample (defining (...) the thermodynamical average) is 

Pj(q) = (6(q- q r)) , (4) 
and averaging over samples (denoting with TTT the average on the disorder) we define 



P(q) = Pj(q) 
The overlap Binder parameter is defined as 



9 =2 



The spatial overlap-overlap correlation function is 



a 
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(5) 



(6) 



(7) 



We will denote by Cj (or C(J)) the q — q correlation at distance j, averaged over different 
site couples whose distance is j. 

It is also possible and interesting to define a correlation function where the overlap 
is fixed to a given value: one selects couples of equilibrium configurations (the two real 
replicas that constitute our system) with a given fixed overlap q and compute the correlation 
function among these configurations. We will denote these correlation functions with the 
symbol C q (x). Obviously 



C(x) = Jdq P(q) C q {x) 



i.e. C(x) is the sum of C q (x) weighted with the static probability distribution of the 
overlaps. Finally we will define the connected correlation functions, by using the symbols 
C q (x) and C, 



C q (x) = C q {x) - q 2 , 
C{x) = C{x)-W). 



(9) 



Some crucial properties of the propagators, as computed in the RSB framework, will be 



useful in the following (see [15 and references therein). At tree level one finds that (in 
dimension D < D" = 6 there are corrections, and the contributions due to the anomalous 
dimension has to be included) 



x 



C q (x) OC 
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if q = 0, 
if < q < q E A 
if q = q EA . 



(10) 
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We have based the analysis of the numerical results of our numerical simulations on finite 
size scaling techniques. When changing the temperature T and the lattice size L the 
overlap susceptibility x = V(q 2 ) scales, in all generality, as 

X {L, T) = Li (h x (L^(T - T c )) + L-"h 2 (l x ' v {T - T c j) + 0(L- 2 ")) , (11) 
and the overlap Binder parameter scales as 

g(L, T) = h (l}' v {T - T c )) + L~ u f 2 (L^{T - T c j) + O^) , (12) 

where u is the exponent that determines the corrections to the scaling (it is the derivative 
of the /3-function at the critical coupling, flnj), and fi, f 2 , hi and h 2 are universal functions. 
In particular this effect is important for the overlap Binder parameter at T = T c , and gives 
scaling violations at the infinite volume critical point. One has that 

g(L,T c ) =/l ( ) + L-/ 2 (0) + O(L- 2 -) . (13) 

This effect explains why there is not an unique crossing point of the Binder cumulant 
curves (see, for instance, reference ||18||), and has to be considered in the analysis of the 



numerical data. A one loop field theoretical computation [£L6] gives 



u = 6 - D , (14) 

i.e. in D = 3 at first order in perturbation theory one expects uj — 3. Of course this number 
will be modified by renormalization, but we expect the correct result not to be far from the 
previous guess. In a related theory, the four dimensional site percolation, also described 
by a generalized 3 theory, an exponent for the corrections to scaling uo ~ 1.13 has been 
determined numerically. That has to be compared to the naive (one loop) expectation, 
uj = 2, and to the three loop value (with Pade resummation), to = 1.52 []19 . In spin 



glasses the computer power available at present does not allow an accurate estimate of 
the correction to scaling exponent. 



3 Numerical Methods 

We will base our discussion on large scale simulations of three dimensional systems of 
linear size L — 4,6, 8, 10, 12 and 16. The (earlier) simulations on smaller lattices use the 



tempering updating scheme ||20|| , while (more recent) simulations on larger lattices use the 
parallel tempering updating scheme [|21|| , which turns out to be very effective (see [E3] for a 
review of improved methods, and [0 for recent work relevant for spin glasses). In all cases 
checking that the systems were, as far as we could establish, fully thermalized, has been 
the first of our worries. 

We have used the APE-100 parallel supercomputer [23]|. Some of the first runs have 



used the tube version (with 128 processors), while the bulk of the simulations have been 
run on the tower version, with 512 processors, with a peak performance of 25 Gflops. On 
the tower our codes has a sustained speed of 6 Gflops, and updates 200 million spins per 
second (even if spins can only take two values, the model is based on Gaussian couplings, 
and floating point arithmetics is needed). 
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The simulated tempering method 0, (together with multicanonical approaches, 
see [0 and references therein) is a quantum leap for simulations of systems with a complex 
phase space. One can now thermalize, in the broken phase, systems that could never be 
studied with normal Monte Carlo. Normal tempering works quite easily up to L = 10 at 



T ~ 0.70T C (even if the fine tuning of the g constants of |2(J is somehow cumbersome). 
For larger lattices, L — 12 and 16, we have found the use of parallel tempering |2lJ| 
mandatory. As a negative caveat we feel like adding that using 2 months of CPU of our 
25 Gflops computer we have not been been able to thermalize a L = 24 lattice in the low 
temperature phase (down to 0.7T C ). 

The common characteristic of tempering and parallel tempering is that the temperature 
becomes a dynamical variable. In the simulated tempering we propose to update the 
temperature of the system at run time, after one or more usual Monte Carlo sweeps of all 
the spins. This fact makes possible to go from the paramagnetic phase to the spin glass 
phase and back; these changes enable the system to avoid the high free energy barriers that 
separate pure states, making possible to escape from the dynamical traps (the metastable 
states). 

The parallel tempering is far easier to program than the simulated tempering. While 
in the simulated tempering we need to estimate with good precision, in a series of thermal- 
ization runs, the relative free energy of two systems at two contiguous temperatures (in 
the set of the allowed T values), in the parallel tempering method this problem is solved 
by the dynamics itself (i.e. the method generates, at run time, the correct free energies). 



For more details see references |2l], y, |22 |. For sake of completeness we give here the 



scheduling of the two approaches. In the simulated tempering method the process is: 

1. We perform iV TERM iterations at a fixed temperature Ta using the Metropolis algo- 
rithm. 

2. We perform AfpE steps with the Metropolis algorithm at the same temperature Ta 
that will be used to compute a first guess for the relative free energies. 

3. We repeat the steps 1 and 2 following an annealing procedure (i.e. we start with 
the highest temperature and finish with the coldest using the final configuration of a 
given temperature as the initial configuration for the next lower temperature). 

4. We perform five cycles of Metropolis plus simulated tempering runs (five series of 
one Metropolis sweep of all the lattice spins at fixed T plus one trial T update), at 
the end of each we tune the values of the relative free energy to make the system 
visiting each allowed T value for the same time. The total number of steps (for the 
five cycles) is A^jfe 

5. Finally we perform Admeasures steps, measuring the interesting physical quantities. 

We report in table (|l|) the values of the parameters used in the simulated tempering 
runs. 

As we already said parallel tempering is simpler. Here there are no free parameters 
to be fine tuned. We first run Aterm thermalization sweeps of the Np copies of the 
system (each at a different T value), using from the start a Monte Carlo sweep for the 
spins followed by trial temperature sweeps among the copies. After thermalizing we run 
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Table 1: Parameters of the simulated tempering runs. 
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Table 2: Parameters of the parallel tempering runs. 



-^measures sweeps, like the ones we just described, during which we measure and average 
the physical quantities. We show the parameters of the parallel tempering in table (0). 
For the L < 12 lattice we have used temperatures ranging from 1.3 to 0.7 with a step of 
0.05, while that for the L = 16 lattice we have taken T ranging from 1.8 down to 0.7 with 
a step of 0.05. 

As we said checking thermalization has been one of the crucial issues of this work. We 
have checked the following facts: 



1. In the both updating schemes the time that a system spends with a given temperature 
must be independent of the temperature. By construction this fact holds only at 
equilibrium. All our runs verify this property. 

2. The acceptance factor for the temperature update has been monitored, and kept in 
the range 0.2 — 0.5. 



3. A set of relations among expectation value has been recently proven by Guerra |24 



see also [25. 26]). These relations are satisfied at equilibrium. All our measurements 



satisfy them with a small error. 

4. The single sample probability distributions of the order parameter Pj(q) have to be 
symmetric for q — > —q. We have verified this symmetry. 

5. We have monitored the growth of the spin glass susceptibility and the overlap Binder 
cumulant. The times that these quantities need to reach their plateau values are 
estimates of thermalization times. We have kept them under control. 



4 The Binder Cumulant and the Susceptibility 

We will start by discussing our analysis of the overlap Binder cumulant, supporting the 
existence of a spin glass phase transition [l], [|. We show our data in figure (0). A magnified 
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view is in figure (§): it makes clear the crossing (the signature of a phase transition) between 
the L = 4 and L = 16 curves and the L = 8 and L = 16 curves. From figure (0) a first 
(qualitative) guess for the critical temperature is 

T c ~ 0.98 ± 0.05 . (15) 

These results support and improve the results of |1[, where the crossing was first ob- 
served. Here we offer evidence for a split of the curves in the low T-region down to T ~ 0.7T C 
(improving over the single point at T ~ 0.9T C of |T|]). We notice that the splitting is much 
smaller in the spin glasses than in usual ferromagnets because also in the infinite volume 
limit, if replica symmetry is broken, the Binder cumulant is not identically equal to one in 
the low temperature phase, but it is given by 

_3 1 JdqP(q)q 4 
9 2 2(/dgP(g)g 2 ) 2 ' 1 ' 

The small dependence of g on L reflects the mild dependence of P{q) on the size. It is 
evident that g does not go to 1 at T < T c (see fig. (|3])) implying that the function P(q) 
does not tend to a delta function in the infinite volume limit in contradiction with the 
droplet model predictions ||. 

The second information that we can obtain from the analysis of the Binder cumulant 
is the value of the v critical exponent. We have estimated the value of the exponent v by 
fitting the derivative of g in the high temperature region with the technique introduced in 



27]. We introduce the function / as the derivative of the Binder cumulant g with respect 



to T at fixed value of g = g : 

An i 

~ aL~ . (17) 



f(9o,L) = 



T :g(T )=g 

We have computed the derivatives of the overlap Binder cumulant in the range from g = 0.6 
to g = 0.8 (the crossing point is close to the value 0.75). In this interval we have fitted the 
Binder cumulant curves (at fixed L) using fourth and fifth order polynomials (in order to 
control the systematic effects due to the fitting procedure). To cross check consistency we 
have also used a derivative both with respect to the temperature and with respect to (3: 
the two procedures are equivalent and they should produce the same result. 

In a large region with T > T c close to the critical point such an estimate turns out to 
be, as hoped and expected, independent of the go value that has been used. We plot in 
figure (|) the values of -(go) as a function of g (obtained fitting the data with a fourth 
degree polynomial both in T and in (3). We expect the fit to break down when going too 
close to g c because of the large error involved. 

The fit to the power behavior of eq. ( P~T|) works well in the range of go going from 
g = 0.6 to g = 0.68 (i.e the % 2 /DF ~ 1, DF stands for degrees of freedom): for g > 0.69 
we get to close to the crossing and the fit becomes unfaithful. This happens both for the 
f3 and the T fit. The error looks bigger in the (3 interpolation. The final value is 

v = 2.00 ±0.15 . (18) 

Let us stress the "locality" of this method. We perform the derivative closer and closer the 
critical point. When we determine a plateau we know we are entering the critical region 
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Figure 1: The overlap Binder cumulant versus T for L = 4,6, 8, 10, 12 and 16 (curves 
from top to bottom on the right part of the figure). 




Figure 2: A magnified view of figure (TO. In the upper window crossing between the L = 8 
and L = 16 lattices. In the lower window crossing between L = 4 and the L = 16 lattices. 
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Figure 3: The overlap Binder cumulant at fixed T = 0.7 versus 1/L. 

where scaling violations are small. In reference [0, Kawashima and Young simulated the 
±1 version of the model using an annealing procedure in order to thermalize the system. 
The larger lattice that they simulated was L = 16. They found v = 2.0 analyzing the 
Binder cumulants curves. Another value of v was obtained analyzing the scaling of P(q) 
near the critical temperature, v = 1.6. (obtained also (3/v). Finally they reported as final 
value v = 1.7(3). 

Instead of study the scaling properties of the overlap probability distribution we will 
follow another way. We will compute, using the Binder cumulant data, the infinite volume 
critical temperature, and using the power law scaling of the susceptibility at the critical 
point we will get the critical exponent j/v. 

The Binder cumulant also gives us quantitative information about the value of the 
critical temperature. One can see that the temperature where the Binder cumulant 
takes a given value, that we call g\ = g(L, T C (L, gi)), close enough to g c , scales as 

T c (L, 9l ) =T c (w) + aL- 1/u . (19) 

Such a value of the infinite volume critical temperature and of the critical exponent will 
depend on g\ if it is not close enough to g\. when taking the infinite volume limit we have 
to keep gi close enough to g c . Inside the critical region we will find that asymptotically 
the critical temperature and the critical exponent become independent of g±. For instance 
the three parameter fit using gi = 0.68 give us 

T c = 0.99 ±0.07 . (20) 

Nevertheless the determination of v has a very large error. In order to reduce the eror bars 
in T c we fix v to the value given in Eq. (|18|): v = 2.00(15) Fixing v = 2 in ([19|) we find 
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Figure 4: ^(<?o) versus g . Triangles for fits over T, squares for fits over (3. 

T c = 0.92(3) for g x = 0.65; T c = 0.95(4) for g x = 0.68 and T c = 0.92(3) for g x = 0.70. For 
gi > 0.70 the quality of the fit is poor. Taking into account all these numerical results and 
the u's error bars [|32| we finally quote the result 



T r = 0.95 ± 0.04 



2.00 ±0.15 . 



(21) 



We end this quantitative study with the standard plot of the overlap Binder cumulant, 
in figure (|5p: we plot g versus the scaling variable L 1 /^(T— T c ), using v = 2.0 and T c = 0.95. 
All the points (from different lattices and temperatures) go to the same curve giving us a 
very good scaling plot. 

After locating with precision the location of the critical point we can study the diver- 
gence of the spin glass susceptibility. The susceptibility diverges with a power law with an 
exponent given by -. Fitting our data at T c = 0.95 we obtain 2 = 2.36. We must take into 
account the error bars on the critical temperature. At T = 0.95 + 0.04 we find 2 = 2.30 
while at T = 0.95 — 0.04 we have 2 = 2.42, so we quote 



7 
v 



2.36 ±0.06 



(22) 



These exponents are in good agreement with those found by Kawashima and Young [p]] for 
the J = ±1 Ising spin glass, v = 1.7±0.3 and £ = 2.35±0.05, and with the £ = 2.37±0.04 
found by Berg and Janke in 0. 
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Figure 5: g versus L~(T — T c ), with u = 2.0 and T c = 0.95. In the precision given by 
the statistical error all the points (from different lattice sizes) collapse on the same curve. 
Triangles for L — 8, squares for L = 12 and pentagons for L = 16. 



5 Analyzing P(q) 



The issue of the properties of the probability distribution of the overlap is controversial. 
The P(q) is a crucial ingredient of the RSB picture: its equilibrium non trivial shape 
implies a whole series of unusual patterns of behaviors. 

In the droplet model || P(q) has very simple properties: it is composed by two Dirac 
deltas at ±gEA, where qea is the Edward-Anderson order parameter. 

Moreover we can cite the approach of reference || about another definition of the 
overlap (the so-called window overlap) which must have a droplet model like probability 
distribution. 

However in reference [f| has been shown that both overlaps, the first one computed in 
all the lattice and the second one, the window overlap, computed in a small region around 
the origin (following the prescriptions of reference U), provide the same picture of the 
probability distribution P(q), supporting both definitions of the overlap the RSB picture. 

We will show below, that P(q) does not follow the shape predicted by the droplet model 
||, following qualitatively the properties predicted by Mean Field. The average properties 
of the function P(q) have been discussed in reference f28| . The results are shown at the 
lowest temperature in figure (Q). The position of the peak (q M ) was fitted as (at the lowest 
temperature we considered, T = 0.7) 



qu 



(0.70 ± 0.02) + (1.6 ± 0.7)L- (1 ' 5±0 - 4) 



(23) 
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Let us discuss how the system develops with increasing L a Dirac delta function at q — 
q EA ~ 0.7. We define 



where q , m ax(-^) is the overlap value where Pidql) is maximum. If in the infinite volume 
limit the system develops a Dirac delta function at qEA = 0.7 A(L) should not depend on 
L (at least for L large enough). The weight of the Dirac delta will be 2 A (if the delta is 
approximated by a symmetric function). 

For instance at T = 0.7 we find A{A) = 0.25 ± 0.03, A(6) = 0.27 ± 0.01, A(8) = 
0.26 ± 0.01, A(W) = 0.22 ± 0.05, A(12) = 0.27 ± 0.03 and finally A(W) = 0.25 ± 0.04. We 
can assume that the previous values are roughly constant (A ~ 0.26 ±0.02), and we expect 
that the weight of the Dirac delta at qEA = 0.7 will be, if the delta is approximated by a 
symmetric function, 0.52 ±0.04. This implies that it is consistent to assume that the peak 
will become a delta function in the infinite volume limit. In any case the weight must be 
substantially larger than 0.26. 

Moreover this result supports the existence of a continuous part of P{q) between 
and qEA = 0.7, that carries the missing probability 0.48 (if the delta is approximated by 
a (almost) symmetric function), in order to make Jq d\q\P(\q\) = 1. This analysis gives 
further support to the existence of a replica symmetry breaking phase transition in the 3D 
EA spin glass. 



Figure 6: P(q) probability distributions at T = 0.7 for all the simulated lattice sizes (right 
to left):4, 6, 8, 10, 12 and 16. 

According to the replica predictions the function Pj(q) changes dramatically from sam- 
ple to sample. A smooth dependence on q is obtained only after average. We show in 
figure ((?]) three different functions Pj(q) for three different samples. One can see that 




(24) 
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Figure 7: Non normalized Pj{q) (histogram) for three different realizations of the quenched 
disorder, L = 16 and T = 0.7. 



different samples have a very different behavior. In the first distribution on the top we 
have a sizable contribution for q ~ 0: these kind of samples contribute to the low q part 
of P{q) (exactly like it happens in mean field: only some samples carry weight at q ~ 0, 
while some samples have a zero probability for q ~ overlaps). The sample in the middle 
has five clear maxima, one in q = and four at high q values. The lower sample has 4 
very sharp maxima (here the similarity with a 5 function starts to be clear), and close to 
no weight at q ~ 0. Notice that the level of symmetry of these plots for q <-» — q is a good 
check of how good the thermalization has been. 

Here we will compare in quantitative way the detailed predictions of the sample to 
sample fluctuations of the function Pj(q) with the very detailed predictions of replica 
theory |3D|, and we will find a remarkable agreement. 



From the Pj(q) of (f|) we define the integrated probability distribution 

xj{q s ) = f S dq Pj{q) = 1 - yj(q s ) , (25) 
J o 

i.e. the total probability that the overlap in a given sample is smaller than q s . In the 



following q s will be kept fixed. The probability distribution of xj(q s ), that following |30 
we denote by Il qs (xj), is also a random variable, and it depends on the disorder realization 
J. In order to get a better statistical signal (for reasons that will become clearer in the 
following) it is useful to define the integrated probability of Il qs (xj) 



n<(r)= / dzU qs (z) . (26) 

We define also 
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x(q s ) = xj(q s ) = [' dq P(q) = 1 - y(q s ) . (27) 
Jo 

With our numerical simulations we compute Pj(q) for every realization of the quenched 
disorder J. We select a fixed value of q s , and for each realization J we compute the 
probability to find at equilibrium an overlap smaller than q s . The answer to the question 
"for which percentage of disorder realizations such a probability is smaller than r ?" is given 
by II< (r). This procedure is repeated for a set of q s values (in the numerical work we have 
used q s = 0.1, 0.2, . . ., 0.9, see below). 



The properties of this distribution have been carefully studied in the past [29, 30 



the framework of the Mean Field RSB theory. An analytic expression can be obtained for 



it; indeed one can prove that |30| 



Av e - fa( „) = l-»(«.)W» (7fe) f (28) 

where U qs (yj) = fl gs (1 — xj) (i.e. it is the probability distribution of the random variable 
yj(q s )), and 



f 1 dz U qs (z)eV^ . (29) 
Jo 



9g 3 W 

The D a are the parabolic cylindric functions, and x(q s ) has been defined in 

The function H qs (z) has many interesting properties. We will concentrate on some of 
the asymptotic properties of H q3 (z) for z close to 1. It has been proven that the previous 
formulae imply that: 



U qs {z) - (1 - z)^- 1 as*->l, (30) 



or equivalently 



fl qs (z) ~ z x{qs) - 1 asz^O. (31) 
Finally, in terms of the integrated probability distribution 

ft<(z) ~ z x{qa) asz-^0. (32) 

The integrated probability distribution from to 2 (for a fixed reference value q s ) goes to 
zero (when z — > 0) as a power law of z with exponent x(q s ). What we are looking for is 
the probability of finding a function Pj(q) which is very small in the region < q < q s (i.e. 
Jq s dqPj(q) < z), and this probability goes to zero as a power of z. This behavior is quite 
peculiar, especially in the region where x(q s ) is small. Indeed we find that fl^(z) remains 
a quantity of order 1 as far as we stay in the region 

z > e - ^) . (33) 
When x(q s ) is small the probability distribution of z becomes very intermittent 



z = x(q sj 

gSIM) = g-jj^j . (34) 
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Where logx means, in the whole paper, the natural logarithm of x. Before using our 
numerical simulations of the 3D model to check if we find a similar power behavior, it is 
convenient to look more closely to the analytic predictions. In the RSB solution of the mean 
field theory one can directly inspect the solution fl28p to exhibit the power law behavior. 
However the computation of the function itself is rather lengthy because it involve the 
evaluations of inverse Laplace transforms. Let us illustrate here a different way to obtain 
the correct answer, which is based on probabilistic techniques and has also the advantage 
of giving a further insight in the features of the RSB solution. 
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Figure 8: log(n<(x)) versus log(x) in the mean field approximation (log means the natural 
logaritm). Here q has been fixed to a value such that x c = 0.2 (see text). The clear power 
law behavior is well fitted by an exponent equal to x c . 



We recall that in the RSB theory the ultrametric structure of the states implies that 
it is possible, by fixing a given value of the overlap q c , to break the set of all the pure 
states in clusters p9| . Two pure states belong to the same cluster if their mutual overlap 
is bigger than q c . This is equivalent to fix a value of x c = x(q c ), since the relation between 
q and x is monotonic (see equation (|27|)). Each cluster is characterized by a weight wj , 
where / is the index of the cluster and J means that we are working in a given disorder 
realization, that is nothing but the sum of the weights of all the pure states belonging to 
the cluster. The probability of finding an overlap greater than q c (the reference value in 
the construction of the clusters) turns out to be 



yj(q c ) 



xj(q c 



w 



J\2 



(35) 



where the index I is running over all the clusters 33 
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A theoretical analysis tells us that the probability distribution of the weights of the 
clusters can be generated by the following algorithm. We consider a system where N 
clusters are allowed (eventually iV will go to infinity). We assign the weights wi to the N 
clusters as follows 

Wl = Ae~ pFl , (36) 

where Fj is the extensive free energy of cluster /, and A is a (F-dependent) normalization 
constant such that J2i=i,n w i = 1- The extensive free energies are mutually independent 
and they are distributed according to 

f Bpf 3x (lc) F I if F < Fat 

P(F I ) = \ Be * 1 ' (37) 

[ otherwise . 

where B is a normalization factor. With these weights we can construct the quantity 
1 — J2i w j an d compute the probability distribution of these sums (the II 9c ) or the integrated 
probability distribution IT5(r). The result for the weights does not depend on Fm, which 
we can take equal to 0. If we use this algorithm in the limit where N — ► oo we obtain the 
correct probability distribution. 

This algorithm can be easily implemented numerically. The probability distribution of 
of the quantity u>i = exp(— (3Fi) is given by 

P Q (UJ) OC -, — t—j (38) 

can be now generated by extracting iV weights w using the random numbers p, uniformly 
distributed in (0,1): w = A-^j^. With these weights we can construct the quantity 

w j an d compute the probability distribution of these sums (the Tl qc ) or the integrated 
probability distribution Il< (r). 

Very good results can be obtained also for moderate values of iV (i.e. iV > 10 if we do 
not consider a value of x close to 1. We show in in figure (El) the result of this algorithm for 
a choice of q c such that x c = 0.20. As it should be the data follow a very accurate power 
law in the range of log(x) (natural log) going to —25 to — 1, with an exponent of 0.2. 

We plot the results for n^(r) from the numerical simulations of the 3D model in figure 
(§). The power law behavior is clear, and consistent with the results of mean field theory. 
The exponent x^r,{q c ) does not coincide exactly, as a function of q c , with the function 
x(q c ) (calculated by integrating directly the equilibrium P(q)): we plot the two functions 
in figure (fLOp . 

We can obtain more insight about P(q) by analyzing the moments of the tl qc (r) prob- 
ability distribution. We consider the relations P9fl , valid in the RSB solution of the mean 
field model, 
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where 



. , 1 . . 2 2 . . 
%2 {q c ) = -xi{q c ) + g2?i(g c ) , (39) 

x 3 (q c ) = (3 + 7m((fc) + 5x 1 (g c ) 2 ) , (40) 
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Figure 9: log IT^ (x) versus logic (log means natural logarithm) from the results of numer- 
ical simulations of the 3D Edwards Anderson spin glass. Here q = 0.3 (which corresponds 
to x c = x(0.3) ~ 0.15), L = 16 and T = 0.8. 



x n (q c ) = J dr r n fl qc (r) . (41) 

We plot in figure (|TT|) ^ versus x\ for the L = 8, 10 and the L = 16 lattices at two different 
temperatures: T = 0.8 and T = 0.7. We also show the straight line prediction of mean 
field theory (|3~9p. We note the numerical data do not depend on the temperature. The 
figure also shows that for x > 0.5 there is a good agreement among the mean field RSB 
solution and the numerical data for the 3D model. For x < 0.5 data on larger lattices are 
becoming closer to the mean field result. Recent analytic and rigorous work p|, 



is 



focusing on relations of the type (|39|), and on the relations among expectation values with 
different number of replicas. 

We have obtained similar results analyzing x^jx\ against x\. Our numerical data are 



plotted in figure ([12]) together with the Mean Field prediction of equation (f40D). 

We have seen that the replica predictions are in good agreement with the numerical 
data. There are however two problems: 



• the relations (^) and ( f4"0"D seems not to work at small X\\ 

• the value of the dynamical exponent £3d(<? c ) is higher that what it should be. 

It is not clear from our data if these discrepancies are a finite volume effects or if they 
disappear when the volume goes to infinity. We want to point out that these discrepancies 
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Figure 10: The exponent x^o{Qc) obtained on a L = 10 lattice (squares) and on a 
L = 16 lattice (triangles) at T = 0.8. The continuous line is the equilibrium function x(q), 
obtained by integrating the equilibrium probability function, P(q), using equilibrium data 
on a L = 12 lattice. 
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Figure 11: ^ versus x\. For L = 10 T = 0.7 (triangles), and T = 0.8 (squares); for L = 8 
T = 0.7 (pentagons), and T = 0.8 (hexagons); for L = 16 at T = 0.7 (heptagons) Triangles 
and pentagons are connected with continuous lines. The straight line is the mean field 
prediction. 
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Figure 12: ^ versus x±. For L = 10 T = 0.7 (triangles), and T = 0.8 (squares). For 
L = 10 T = 0.7 (pentagons). The upper continuous line is the mean field prediction. 



are just what we would expect from finite size effects. Indeed in a finite volume we could 
expect that 

P ((l) = S W a W pf(<l ~ 9a,/?» A «,/3) » ( 42 ) 

where the function f(z,A) is essentially different from zero for \z\ < A. In the limit 
A — > we should have F(z,A) — > It is obvious that delta function only appears 

in the infinite volume limit and the quantity A is likely to go to zero in the average as 
a power of L. In this situation the quantity Xj is always of order q when q goes to zero 
(the function Pj(q) being bounded) and therefore the expectation of x 2 j is of order q 2 . 
Consequently in the region where q < A we must have that the ratio — must go to zero 
with x. It is remarkable that this behavior sets in only at small values of X\. 

The same argument can be applied to the tails of the probability distribution. For sake 
of simplicity we consider the case where 

/(*,AH*^> (43, 

It is evident that the the value of P(q) will be affected by the presence of states having 
an average overlap different at most A from q. A simple computation show that in the 
approximation where A does not depend on the states a and f3, we have that 

x 3D {q c ) =x{q- A), (44) 
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which is roughly what we observe, with A « 0.15. This value is of the same order of 
magnitude of the width of the function P(q) at the peak and also is of the same order 
of magnitude of the region, in x±, in which the relations (|39|) and ( f40|) do not work. A 
much more accurate analysis is needed to decide if these small discrepancies arise from the 
mechanism outlined here and consequently disappear at large volumes. 

6 A Simple Function 

Let us consider the function 



0(1 -<M», (45) 

that we will call F. In the RSB solution of the mean field model one can compute its 
behavior. Usual arguments imply that F = x, X being the magnetic susceptibility. In the 
high temperature phase P(q) is a delta function, and (\q\) = 0: everybody agrees that this 
result also holds in finite dimensional, realistic spin glasses. 
In the SK model one finds that 



T/T c 



(46) 



so that F is exactly constant. 
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Figure 13: F(T) versus T for L = 4 ... 16. 



This behavior will be in any cases changed by renormalization at finite D; (\q\) does 
not go to zero but it is proportional to (1 — T/Tc) 13 , , where the the exponent (3 is equal to 
|z/(l— 77), which in three dimensions is not far from 1 (our best estimate is (3 = 1.30±0.12). 
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It is interesting to study how much in 3D this function stays close to a constant in the 
cold phase also because it should equal to the magnetic susceptibility. We plot F in figure 
(|I~3|) for different values of L: the function turns out to be remarkably constant for T < T c , 
where q has a highly non-trivial behavior. 



7 The Correlation Functions 

In reference || the q = ergodic component of the overlap correlation function was com- 
puted by using an off-equilibrium dynamics. The off-equilibrium approach was used since 
it is difficult to thermalize large lattices in the cold phase, and an extrapolation procedure 
was used to deduce the equilibrium decay rate on the relevant distance scales. The power 
law decay of this correlation function gives strong support to the validity of the mean field 
picture in finite dimensional systems. Here we are able to compute the same correlation 
functions directly at equilibrium, and to show that we find a completely consistent result, 
giving in this way even larger support to the validity of the RSB Ansatz for describing 
realistic spin glasses. 

We recall the definition of the equilibrium non-connected correlation function. For each 
pairs of configuration o and r we define the correlation C(x) as 

c ( x ) = Ts E a ( x + yMvM* + y) r (y) • ( 47 ) 

L y 

It may be convenient to recall that C(x) does not change under a global flip of the a or r 
configurations. If we restrict the statistic over configurations that have a given overlap q, 
the correlation function that we will obtain will be denoted by C q (x) or C(x\q). 

The upper curve in figure ( |Hp is the equilibrium correlation function at T = 0.7 on the 
L = 16 lattice. The two lower curves are off-equilibrium q = correlation functions from 
0: in the lower curve the system was suddenly quenched from T = oo to T = 0.7, while in 
the second one a slow annealing procedure was used to bring T down to 0.7. In both cases 
the lattice size was 64. The third curve from the bottom has been computed, by the same 
configurations used in the upper curve (equilibrium runs at T = 0.7 and L = 16), including 
only couples of configurations where q < 0.01. This is the real equilibrium q ~ correlation 
function, where no extrapolations were needed. The agreement with the result of the off- 
equilibrium runs is very good: only at x ~ 8, close to the center of the periodic lattice, we 
see an expected discrepancy (the lattice of the off-equilibrium runs had a linear size of 64, 
and only close to L = 32 boundary conditions became relevant). All the three curves are 
well compatible with a power law decay of the form x~ x with A ~ 0.5 in agreement with 



the predictions of |L5] . 



As we have seen the propagator at fixed q is equal to q 2 plus a a term which goes to 
zero at infinity (see eq. fllPP). Therefore at large distances x 

Pc{C{x)) = P(C(X)2 1 ) , (48) 
(2C(x))5 

producing a double peak structure, one peak corresponds to the delta function at qteA and 
the other peak comes from the singular Jacobian (of the transformation C oc q 2 ) at zero 
overlap. 
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Figure 14: The lower curve is the infinite time extrapolation of the non-equilibrium 
correlation function Cd(x\q = 0) obtained by a sudden quench (L = 64). The second curve 
from the bottom is Cd(x\q = 0) obtained by a slow annealing (L = 64). The third curve 
is the equilibrium correlation function computed with the constraint \q\ < 0.01 (L = 16). 
The upper curve is the full equilibrium correlation function, including all configurations (L 
= 16). 
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If we naively assume that C q (x) = a(x) + b(x)q 2 both peaks are present also for finite x. 
If we consider a more complicated dependence of C q (x) at finite x the leftmost peak may 
disappear and one could be left with a non divergent structure. A divergent peak should 
however be present for large x. 

In order to justify the equality among the dynamical result and the static results at 
q = we can intuitively argue as follows. We define an "effective" potential [0 given 
by — log P(C(x)). In an off-equilibrium simulation with random initial conditions all the 
components of the propagator start at zero. Using the effective potential — log P(C(x)) 
it is clear that the initial C(x) — values drifts to the closest minimum, that is given by 

Cmm(x), i-G 



lim C d (x,t\q = 0) = C min (x) 



(49) 



where Cd(x,t\q = 0) is the dynamic correlation function computed at time t and distance 
x using a sudden quench from infinite temperature down to T < T c (so that the initial 
overlap is zero, fact that we denote with q = in the formula). 
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Figure 15: Histogram of the equilibrium correlation function at distance x — 1. We have 
marked with a vertical line the value of Cd(l, t\q = 0) extrapolated to infinite time [Of. 



The equilibrium probability distribution of C(x) should show a double peak structure. 
The smaller peak should be located at C m i n (x), and this value must be equal to the off- 
equilibrium value Cd(x,oo\q = 0) in the infinite time limit (i.e. the lower curve of figure 
([14])). We show these histograms in figures ([15]) and (16) for x = 1 and x = 3 respectively. 
We have also marked with a vertical line the value obtained from the off-equilibrium sim- 
ulations. 
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Figure 16: Histogram of the equilibrium correlation function at distance x = 3. We have 
marked with a vertical line the value of Cd(3,t\q = 0) extrapolated to infinite time ||. 



Qualitatively we see that both histograms have a clear peak (for large values of C, 
that correspond to C max ) and a second maximum or flex point (for lower values of C, that 
correspond with C m i n ). We have also shown the value of the infinite time extrapolation of 
Cd(x,t\q = 0) (the vertical line in the plots). The same pattern holds up to x = 5 (that is 
the larger distance that we analyzed in the dynamical runs of O). It is clear from figures 
pp|) that the value of C(x\q = 0) is located close to the second maximum or flex 



(15|) and 
point. 

This fact strongly supports the correctness of the off-equilibrium approach for the com- 
putation of the q = component of the propagator (including the effectiveness of the 
extrapolation procedure) and confirms that it is possible to compute equilibrium expec- 
tation values in off-equilibrium simulations. Moreover the double peak structure of the 
correlation function probability distribution survives in the infinite volume limit since the 
dynamical expectation values have been computed on a very large lattice (infinite to all 
practical effects). It is remarkable that the distribution probability of C(x) becomes wider 
with increasing x. 

All these features of the correlation functions provide a further evidence of the existence 
of a non trivial probability distribution of the overlaps, and of the fact that finite size effects 
are well under control. 

Finally we can check our previous estimate of the r\ exponent (j/u = 2 — T], and so 
r] = —0.36(6)). We recall that at the critical point we have 



C(x) oc 



x 



D-2+r] 



(50) 
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This formula always holds at the critical point. 

We show in figure (|T7p the overlap-overlap correlation function near the critical point 
(T = 1.0) for L = 16. Taking into account the points in the interval [1, 5] we have found a 
very good power law fit with exponent 1 + r\ — 0.710(5), giving 77 ~ —0.29 according with 
our previous estimate. Obviously for x > 5 the propagator is seeing the periodic boundary 
conditions and the propagator misses the pure power law behavior. We can remark that 
this early power law behavior that we have found (for x > 1) have been seen, for instance, 
in the four and six dimensional Ising spin glass jl0[ PH . 
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Figure 17: The correlation function near the critical point (T = 1.0) for a L = 16 lattice 
in a double (natural) logarithmic scale. We have marked with a straight line the results 
of a power law fit (taking the points in the interval [1, 5]). The slope is in good agreement 
with the results found by fitting the susceptibility as a function of the lattice size. See the 
text for more details. 



8 Conclusions 

The numerical simulations we have discussed here have allowed to establish some important 
results, and to make clear some relevant issues. Maybe the main point is that the observed 
equilibrium behavior of the 3D spin glass shares the crucial features of RSB solution of 
the mean field theory. We are able to establish that on large lattices, deep in the broken 
phase. 

The study of the overlap Binder cumulant and of the spin glass susceptibility allow a 
precise determination of the critical exponents. The analysis of sample to sample variations 
P(q) allows to exhibit clear power law behaviors in the probability distributions that are 
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typical of the mean field theory (and would have no reason to appear in any theory of an 
usual ferromagnet). Correlation functions are now under control even at equilibrium on 
large lattices, for the whole set of stable configurations and for the zero overlap part of the 
phase space. 

We believe that a very clear picture is emerging. 
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[33] It is easy to obtain this formula. In terms of the weights of the pure states w J a the 
function Pj(q) is given by 

Pj(l) = j:w J a w J p 8( q - q f), (51) 
where a and (3 run over all the pure states. Now, 
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Vj(q c ) = f dq'PAq') = £ w J a w J ,9(qf - q c ) = £ £ , (52) 

qc a. (3 1 ot{I),/3{I) 

where the index / runs over the clusters (defined with the reference value q c ) and «(/) 
and f3(I) run over all the pure states belonging to the cluster I. Taking into account 
the fact that wj = J2 a (i) w a(i) we finally obtain 

vAic) = EK) 2 • (53) 

[34] The partition function is 

Z = J d[0] exp [-nm . (54) 

For a given distance x the value of the correlation function at distance x is a function 
of the field 0: C x = G(4>) = V^ 1 J dycj)(x + y)4>(y), V = J dy is the volume of the 
system. Introducing the factor 1 = / dC x 5[C X — G(4>)} in (|54]) we obtain 

Z = JdC x J d[<f>] 6[C X - G[4>}} exp [-H[<f>]\ . (55) 

We define 

P[C X ] = J d[0] S[C X - G[<f>}} exp [-n[<f>]] (56) 

the probability of finding a given value of C x . This allows to define the effective free 
energy (or effective potential) F[C X ] 

F[C X ] ee- log P[C X ] , (57) 

such that 

Z = J dC x exp[-F[C x ]] . (58) 

In the thermodynamical limit the partition function will be dominated by the minima 
of F[C X ], that justifies calling it an "effective potential". 
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